Skip to main content

Full text of "mit :: ai :: aim :: AIM-295"

See other formats



Artificial Intelligence October 1973 

demo. No. 235 

Berthold K. P. Horn 

The intensity at a point in ^n image is the product of the 
reflectance at the corresponding object point and the intensity of 
illumination at that point. Lie are able to perceive llghtneee, a 
quantity closely correlated ulth reflectance. Hon then do Me eliminate 
the component due to illumination from the Image on our retina? The tuo 
components of image Intensity differ in their epatial distribution* A 
method Is presented here which takes advantage of thle to compute 
lightness from image intensity in a layered, parallel fashion. 

The method is developed for a restricted class of images first used 
by Land in presenting his retinex theory of color. In thle theory the 
problem of color perception is reduced to one of judging black and white 
lightness on three images taken in different parts of ths visual 
spectrum. The method described h*re fills the need for a lightness 
judging procees. 

The theory has implications for potential special purpose hardware 
In Image sensing devices. It should also bs of interest to cognitive 
peychologiete elnce it can explain certain effects observed in the human 
vleu*l system ae well ae predict net* ones. Further, the theory provides 
neuro-phyelologlsts with suggsstlons about the function of certain 
etructuree In the primate retina. 

This report describes research done at the Artificial Intelligence 
Laboratory of the Massachusetts Institute of Technology. Support 
for the laboratory's artificial intelligence research is provided 
in part by the Advanced Research Projects Agency of the Department 
of Defense under Office of Naval Research contract N000U-70-A- 

Reproduction of thle document. In whole or in part, la permitted for any 
purpose of the United States Government* 

1. Revieu. 3 

1.1 Theories of Color Perception. 3 

1.1.1 Tri-Stlmulue Theory. * 

1.1.2 Color Conversion. 4 

1.2 Land's Retinex Theory. S 

1.2.1 Lightness Judging. 6 

1.2.2 fllni-HorJd of Mondriane. 6 

1.3 Uhy Study the One-dimensional Case? 7 
1.3.1 Notation. 8 

1.4 One-Dimensional Method - Continoue Case. 9 

1.4.1 One-Dlmenslonal Continoue Method: Details. 9 

1.4.2 Normalization. 12 

1.5 One-Dimensional Method - Diecrete Case. 13 

1.5.1 One-Dimensional Di ecrete flethod: Details. 13 

1.5.2 Selecting the Threshold. IS 

1.5.3 Accuracy of the Reconstruction. 17 

1.5.4 Generalizations. 17 

1.5.5 Physical Models of the One-Dimensional Discrete Process* 18 

2. Lightness in Two Dimensional Images. 19 

2.1 Tuo-Dineneional Method - Continoue Case. 19 

2.1.1 Applying the Laplaclan to a Mondrian. 19 

2.1.2 Inverse of the Laplacian Operator, 29 
2*1.3 Uhy one can use the Convolutions! Inverse. 23 

2.1.4 Normalization. 23 

2.1.5 Tuo-D i mens I ona I Continous Methodr Detal Is- ♦ - * 24 

2.2 Tuo-Dimeneional Method - Discrete Caee. 25 

2.2.1 Teeselatlon of the Image Plane. 25 

2.2.2 Diecrete Analogue of the Laplacian. 26 

2.2.3 Inverse of the Diecrete Operator. 27 

2.2.4 Computational Effort and Simplification. 29 
2,2*5 Tuo-Dimeneional Discrete Methods Oetatle. 38 

2.2.6 Simplicity of the Inverse. 32 

2.2.7 Iterative Methods of Solution, « 32 

2.2.8 Convergence of Iterative and Feed-back Schemes. 33 

2.2.9 Setting the Threehold. 34 
2.2.18 Some Notes on This Msthod. 35 
2.2*11 Dynamic Range Reduction. 36 
2.2.12 A Frequency Domain Interpretation. 36 

2.3 Physical Mode J e. 37 

2.3.1 A Diecrete Physical Modsl. 38 

2.3.2 A Feed-back Schsme for the Inverse. 48 

2.4 Limitations of the Simple Scheme Presented. 42 

2.5 Computer Simulation of the Discrete Method. 44 
2.5.1 Form of Inveree ueed in the Computer Simulation. 4S 

3. Impl icatione and Concluelone* 47 

3.1 Parallel Image Proceeeing Hardware. 47 

3.2 Cognitive Peychology, *7 

3.3 Neuro-phyeiology. 4 * 

3.4 Conclusion. ** 



I wish to thank Dr. David flarr for a numbsr of most Interesting 
dlacueeione and for encouraging na to commit this msthod to papar. Ha 
first observed that the primate retina appears to have Juet the right 
detailed etructure to implement the functions here described. 

LIGHTNESS* Definition 

The relative degree to which an object reflecte light. 

The Random House Dictionary 

The attribute of object colors by which the object appears to reflect or 
transmit more or less of the incident light* 

Uabster's Seventh Nsu Collaglata Qictionary 


Part 1 Is a review of the relevant Information relating to color 
vleion and lightness. Thie Includes a discussion of the Land ratine* 
model in a form suitable for the developments of the next part. 

In part 2, Land's one-dimensional operation will be extended to 
two-dimensional Images. The method depends on a layered, parallel 
computation suggestive of both biological and artificial 

In part 3 some of the implications are explored and the information 
on the new image processing technique is summarized. 


1. Review, 

1.1 Theories of Color Perception. 

Thar* has always been great interest in how we perceive colors and 
numerous explanations have been forwarded (Newton 1704, Goethe) 1818, 
Young 1820, Maxwell 1856, Helisholt* 1867, Hering 1875). The huaan 
perceptual apparatus Is remarkably eucceeful in coping with large 
variations in the i I lueinatlon. The colors us perceive are closely 
correlated ulth the surface colours of the objecte viewed, deepite large 
temporal and epatial differences in color and intensity of the Incident 
light. Thie is eurprieing eince we cannot eenee reflectance directly. 
The light inteneitg at a point in the iaage le the product of the 
reflectance at the corresponding object point and the intensity of 
illumination at that point - aside from a constant factor that depende 
on the optical arrangment. There must then be soee difference between 
theee two components of Image intensity which allows ue to dlecount the 
effect of one. The two componente differ In their spatial distribution. 
Incident light inteneitg will usually vary smoothly, with no 
di econt inui t ies, uhlls reflectance will have eharp diecontinul ties at 
edges where objecte adjoin. The rsflsctance being relatively constant 
between euch edgea. 


1*1.1 Trl-Stlmului Theory* 

Some facto about how we see color are fairly uell established. It 
appears that us have three kinds of sensors operating In bright 
1 1 lull n*tl on, Hith peak aanai tivi tlee in different parts of the visible 
spectrum. This is why it takes exactly three colors in additive 
mixture to natch an unknoun color. I4h1 le It la a bit tricky to Measure 
the sensitivity curvss of the three sensors directly, a linear transform 
of these curves has been known accurately for some tine (Brlndley 19GB). 
Theee curves, ceiled the standard observer curves, are sufficient to 
allow one to predict color matches made by subjecte with normal colour 
vieion (Hardy 193G1 . 

The elmplest theory of color perception then amounte to locally 
comparing the outputs of three such sensors and assigning colour on thle 
basis (Young 1820, Helmholtz 18671. This however totally falls to 
explain the observed color contancy. Perceived color does not depend 
directly on the relative amounts of light msasured by the three sensors 
(Land 1959, Lettvin 1967). 

1.1.2 Color Conversion. 

A number of attempts have been mads to patch up this theory under 
the rubrics of "discounting of the lllumlnanf, "contraet effect 
adjuetment" and "adaptation". The mors complicated theories are based 
on models with large numbers of paramstsrs which era adjusted according 


to empirical data IHelson 1938 & 1943, Judd 1348 & 1952, Richard* 1971). 
Theee theories are at leaet partially effective in predicting human 
color perception when applied to simple arrangements of stimuli similar 
to those used in determining the parameters. 

The parameters depend strongly on the data and alight experimental 
variations will produce large fluctuations in them. This Is a phenomena 
familiar to numerical analysts fitting curves to data when the number of 
parameters is large. These theories are lacking in parsimony and 
convincing physiological counter-parts. Lettvln has demonstrated the 
hopeleesneae of trying to find flxsd transformations from locally 
compared output of sensors to perceived color {Lettvln 1967). 

1.2 Land*s Retinex Theory. 

Another theory of color perception is embodied In Land'e retinex 
model (Land 1953. 1964 & 19711. Land proposes that the three eete of 

seneore are not connected locally, but instead are treated ae If they 
represent points on three separate images. Processing is performed on 
each euch image separately to remove the component of Intensity due to 

illumination gradient. Such processing is not msrely an added frill but 

ie indi sponsible to color perception in the face of the variability of 

the i I lumination. 


1*2.1 Lightness Judging* 

In essence a judge of lightness processes each image* Lightness is 
the perceptual quantity closely correlated with surface reflectance. 
Only after this process can the three images be compared to reliably 
determir>e colors locally. It remains to mechanize this process. 

It Mould appeal to Intuition If this process could be carried out 
in a parallel fashion that doss not depend on previous knowledge of the 
scene viewed* This is because colors are so immediate, and seldom 
depend on one's interpretation of the scene. Colore will be seen even 
when the picture makes no ssnse in terns of previous experience. Also, 
color Is eeen at every point In an Image. 

1*2.2 lllni-world of tlondrians. 

In developing and explaining his theory Land needed to postpone 
dealing ulth the full complexity of arbitrary scenes* He eelected a 
particular claes of objscts as Inputs, modelled after the paintings of 
the turn-of-the-century Dutch artist Pieter Cornells (londrian. These 
scenes are flat areas divided into sub-regions of uniform matte color* 
Problems such as those occasionsd by shadows and specular reflection are 
avoided In this way. One also avoids shading; that is, the variation In 
reflectance with the orientation of the surface In respect to the sensor 
and the light-source. For Hondrians* lightness is considered to be a 
function of reflectance* 


llondrlans ara usually made of polygonal regione mth straight aides 
— for the development here however the edges may be curved. In the 
world of Mondrians one finds that the reflectance has sharp 
discontinuities wherever regions meet, being constant inaida each 
region. The illumination, on tha other hand, varies anoothlg over the 

1*3 Uhy Study the One-dimensional Case? 

Imagee MrB tuo-dlmeneional and usually sampled at discrete points. 
For hlatoric reasons and intuitive simplicity tha results will flret be 
developed in one dimension, that is with functions of one variable* 
Similarly, contlnoue functlone uill be used at flret eince they aJIOM a 
cleaner separation of the two componsnts of image intensity and 
iiluetrate nore clearly the concepts involved. 

Use Hill be Bade of analogies between the one-dimensional and two- 
dimensional casss as well ae the continous and discrete ones. The final 
procaee discussed for processing image intensities Is tuo-dlmenelonal 
and dlecrete. A number of physical lmplenentatlons for this scheme ar9 

The proceae hIII be looked at from a number of points of view: 
partial dl f ferent ial equat ione, I inear systems, four ter transforms and 
convolutions, difference equations, iterative solutions, feed-back 
schemes and physical models. 


1*3.1 Notation. 

The following notation mi I I be used. 

a* Intensity of incident illumination at a point on the object* 

r' Reflectance at a point of the object* 

p* Intensity at an image point* Product of •' and r'. 

a t r. pi Logarithms of e\ r* and p' respectively. 

d Result of applying forward or differencing operator to p. 

t Result of applying threshold operator to d. 

I Result of applying inverse or summing operator to t* 

Simple derivative operator in one dimension* 

T Contlnous threshold operator, discards finite part. 

1 Simple integration operator in one dimension, 

L Laplacian operator - eum of second partial derivatives. 

G' Inveree of the Laplacian, convolution uith (1/2 t ) log c d/r). 

D*, T*, I*, L* and Gfri Dlacrete analogues of 0, T, I, L and G. 

The output It will not be called lightness since there le probably 
not yet a generally acceptable definition of this term. It is however 
intended to be monotonical ly related to lightness. Note that I le 


related to the logarithm of reflectance, uhile the perceptual quantity 
le perhape more closely related to the square-root of reflectance. 

1.4 Qne-Oimensional Plethod - Continous Case. 

Land invented a simple method for separating the Image components 
in one dimension* First one takes logarithms to convert the product 

Into a sum* This Is followed by differentiation. The derivative uill 
be the sum of the derivatives of the two components. The edges uill 
produce sharp pulses of arsa proportional to the Intensltg eteps between 
regione - uhi le the spatial variation of II I urn (nation uill produce only 
f i nf to values overyuhere. Noh if ono diacarda all finite yatuae, one It 

left ulth the pulses and hsnce the derivative of lightness* Finally one 
undoee the differentiation by simple integration. 

1.4.1 One-Dimensional Continous Method! Details. 

Ue have the following. Let r'UI be the reflectance of the object 
at the point corresponding to the image point x. Let eMx) be the 
intensity at this objact point. Let p* (x) be their product, that ie, 
the intensity recorded In the imags at point x. Note that e' (x) and 
r* (k) %r9 poei tive, 

p* <x) - e' (x) * rMx) 








■ H 

FIGURE 1: Processing steps in the one-dimensional continous case* 


Now let p(x) be the logarithm of p' (x) and bo oni 

ptx) - six) + r(x) 

Note that six) is continous and that r(x) has some finite 

diacont inui ties. Let repcssent differentiation with reapect to x* 

dfx) - D(p(x>) - 0(a(x)) + D(p(x)) 

Now, D(s(x)) Hill be finite everywhere, while 0(r(x)) Mill be zero aalde 
from a number of pulses - which carry ail the information. Each pulae 
Mill correspond to an edge between regions and have area proportional to 
the intensity atep. If now one ■thresholds'* and discards all finite 
parte, one getei 

Km) - T(D(p(x)>) ■ D(r(»)) 

To obtain r(x) one only has to invert the differentiation, that is, 
integrate. Let I represent integration with respect to x f then (I) - 
D andi 

IlKl - I(T(D(p(x))) - rtx) + c 

One can give a convolutions! Interpretation to the above, since 
differentiation corresponds to convolution with a pulse-pair, one 


negative and one positive, each of unit area. Integration corresponds 
to convolution with the unit step function. 

1.4.2 Normalization* 

The reeult is not unique because of the constant introduced by the 
integration. The ztro (spatial) frequency term has been loet in the 
differentiation, so cannot be reconstructed. Thle le related to the 
fact that one does not knoM the overall level of Illumination and hence 
cannot tell whether en object appears dark because it is grey or because 
the level of illumination is loti. 

One can normalize the result If one assumes that there are no light 
•oyrces in the field of view and no floureecent colore or specular 
reflections, Thle is certainly the case for the flondrians, Perhape 
the beet way of normalising the result Is to simply assume that the 
highest value of lightness corresponds to white* or total reflectance in 
the Lambertian sense. This normalization will lead one astray if the 
image doee not contain a region corresponding to a white patch in the 
ecene, but this is the beet one can do. Other normalization techniques 
might involve adjusting weighted local averages, but thle would then no 
longer amount to reconstruction of reflectance. 


l.S Qne-Oimenaional Hethod - Discrete Case. 

So far ue have assumed that tha image intensity uan a continoue 
function. In retinas found In animal a or artificial onss constructed 
out of dlacrata componsnte. iniagaa ara only saitplsd at discrete points'. 
So ona haa to find diacrats analogues for the operations we have been 
using. Perhaps tha simplest ar* first differences and summation as 
analogues of differentiation and integration respectively. This is not 
to say that other approxi mat ions could not be used equally wall. 

To use the new operators one goes through essentially tha sane 
process as before, except that now all values in the differenced imaga 
ara finite. This has the effect of forcing one to choose a threshold 
for the thresholding function. Both components of image intensity 
produce finite values after tha differencing operation. The component 
due to the edges in the reflectance is hopefully quits large compared to 
that due to illumination gradient. One has to find a level that Mill 
suppress the illumination gradient inside regions, uhi la permitting the 
effecte due to edges to reaain, 

1*8.1 One-Dimensional Discrete flethod: Details. 

Let rl be the reflectance of the object at the point corresponding 
to the image point 1. Let s\ be the incident light intensity at this 
object point. Let p\ be their product* that la tha lntenelty In tha 
image at point I. 

- .* 





FIGURE 2: Processing steps In the one-dimensional discrete case, 


Now 1st pj b* th* log of p- and bo on. Lot Oft and Ift 0* !"• operators 
corresponding to taking first differences and summation respect Ively. 
Note that (Ift)"' - D*. 

Pi- •l+n 

d i- p*h- p; (d - d * ( p ,j 
t ; - d; If !d,[ < «. tlio 


ll- > t K - 


1*5.2 Selecting the Threehotd. 

Uhat determines the threshold? It mutt be smaller than the 
emalleet inteneitg etep between regions. It ituet on the other hand be 
larger than values produced by first differencing the maximum 
Illumination gradients. Real images are noisy and the threshold ehould 
be large enough to eliminate this noise inside regions* 

The spacing of the seneor cells must also be taken Into account* 


As this spacing becossv smaller, the contribution due to Illumination 
gradients decreases, uhlls the coaponsnt due to the edges remains 
constant. A limit is reached uhen the component due to illumination 
gradients falls below that due to noise or when the optical properties 
of the Imaging system begin to have a deliterlous effect. In all 
imaging systems an edge Is spread over a finite distance due to 
diffraction and uncorrected abberatione. The spacing of seneore should 
not be much smaller than this distance to avoid reducing the component 
due to edges In the differenced iaags. 

Let u be the radius of the point-spread-function of the optical 
system and h the spacing of ths ssnsor celle. Let g* be the smallest 
atop in the logarithm of raflectanca In tha acana. Than define tha 
affective minimus step asi 

g - g* a mind, h/2u) 

Let a be the largest slope due to illumination gradient and * the root- 
mean-equare noise-amplitude. The noise will exceed a value 3 only «3% 
of the time. Choose the threshold e as follows* 

a < g 

• > a * h 

e > 3\/2 a 


1.6.3 Accuracy of the Reconstruction. 

In thA contlnoue case one can exactly reconstruct tha reflectance, 
aeide from a constant. Us are not eo fortunate here, even If ue Mlect 
a threshold according to the above criteria. This la because the valuee 
at the edges contain email contributions due to Illumination gradient 
and nolee. A alight Inaccuracy In the reconetructlon ulll reeult. Thle 
error le minimized by making the esneor cell spacing very fine, 
optimally of a elze comaeneurate with the optical reeolutlon of the 
device. The effect of nolee can also be minimized by integrating over 

Note that the reconstruction le aore accurate when there are fen 
edgee, elnce It le at the edges that the error effects appear. Ulth 
many edgee the Illumination gradient begins to "show through - . 

1*5.4 General izat lone. 

So far we have dealt ulth constant sensor epaclng. Clearly aa long 
ae the aame epaclng le ueed for both the differencing and the summing, 
the cell epaclng can be arbitrary and has little effect on the 

reconetruction elnce It does not enter Into the equatlone. 


Similarly ue have chosen first differences as the dlecrete analogue 
for differentiation. Ue could have chosen some other uelghted 
difference and developed • eultable Inverse for It. Thle Inveree of 
course uould no longer be eummatlon but cin be readily obtained using 


techniques developed for dealing with difference equetlone (Rlchtmeyer 
1957, Garabedlen 19G4i. 

1*5.5 Physical rtodels of the Ons-Oimensional Discrete Proceee. 

One can invent a number of physical modele of the above operatlone. 
A simple reeistive netuork ulll do for the summation proceee for 
example. Land has implemented a small circular "retina" ulth about 1G 
eensora. This model employs electronic componente to perform the 
operations of taking logarithms, differencing, thresholding and summing. 

Land hee tried to extend his one-dimensional method to Images* by 
covering the image ulth paths produced by a random ualk procedure and 
applying methods like the above to each of these paths. Uhllo this 
produces reeulte. It esens unsatisfactory from the point of view of 
euggeetlng possible neuro-physlologlcal structureet neither does It 
lend itself to efficient Implementation. 

Methods depending on non-linear processing of the gradient along 
paths in the image fall to smoothly generalize to two dimensions, and 
cannot predict the appsarance of images in uhich different pethe reeuit 
In different lightneseee. 


2.' Lightness in Two Dimensional [mages. 

2.1 Tuo-Oi mansions! flathod - Continous CasSt 

Ue naed to extend our Ideas to tuo dimensions In order to deal wl th 
actual images, Thsra are a number of Mays of arriving at tha process to 
be described here, we shall follow tha simplest (Horn 19681. Urn need to 
find two-dimensional analogues to differentiation and Integration. The 
first partial derlvatlvee %r% directional and thus unsuitable alnce they 
Mill for example completely eliminate evidence of edges running In a 
direction parallel to their direction of differentiation. Exploring the 
partial derivatives and their linear combinations one finds that the 
Laplacian operator la tha louest order combination that la leotroplc, or 
rotational ly symmetric. Tha Laplacian operator is of course the sum of 
the second partial derivatives. 

2.1.1 Applying the Laplacian to a flondrian. 

Before investigating the invert ibi I i ty of thie operator* let ue eee 
yhat happene when one appliee It to the image of a Mondrlan. Inelde any 
region one Mill obtain a finite value due to the variation In 
Illumination Intenelty. At each edge one Hi 1 1 get a pulae pair, one 
positive and one negative. The area of each pulse ulll be equal to the 
Intenelty step. 

. -* 




2 j 

FIGURE 3: Applying the laplacian operator to the image of a 
Hondrian figure. 


Thia can beat be esen by considsring tha first derivative of a 
step, nanely a single pulee. If thie ia differentiated again one 
obtains a doubled pulse as described. Since this pulse Mill extend 
along the edge, one mag think of it as a pulse-wail. So each edge 
separating regions m1 1 1 produce a doubled pulse nail. It la clear that 
one can once again separata the coaponant due to reflectance and 
i I I umi nation simply by discarding all finite parts. 

2.1*2 Inverse of the Laplacian Opsrator. 

To co«p)ete the task at hand one than has to find a process for 
undoing the effect of applying the Laplacian. Again there are a number 
of approaches to this problem, we will use the ahortest IHorn 1968). In 
essence one has to solve for p(x,y) in a partial differential equation 
of the formi 

L(p(x,y)) - d(x,y) 

Thie ia Poiason'e equation and It la usually solved inside a bounded 
region using Green's function (Garabsdien 13541 1 

«K,y) - MC{L^« t g)td^^)d{di, 


The form of Gnm'i function G» dependa on the ehape of the region 
boundary Nom If the retina ie infinite all polnta are treated 
similarly and Green* a function dapenda only on tuo para we tare, (J; - x) 
and (*) - y) . Thla poeit tonal indepandanca lapllee that the above 
integral almply beconea a convolution. It can ba ehown that Green'e 
function for thi a caaa lai 

G(£ t *J * *,y) - Q/2 * ) log f (l/r) 

Uhere r x - <£,- x) + (*|- y) 

So p(x 

,y) - / / (1/2") log ( (l/r) * dt^,*| ) d^ d^ 

Thua the Invaraa of the Laplacian operatora la almply convolution ulth 
(1/2*) log <l/r). To ba preciee one haet 

(|-=+ i*£ // d/2"> log c (l/r) * d(^,T ) d^dij - d(x,y) 

Thie ie the two-dimeneional analogue oft 



fit) dt - f(x) 

2*1.3 Why one can use the Convolutional Inverse. 

If the retina le coneidsred infinite one can express the Inverse ae 
a elmple convolution. If the retina le finite on the ether hand one hae 
to uee the acre complicated Creen*8 function foraulation. 

Now conelder e scene on a uniform background whose Image le totally 
contained on the retina. The result of applying the forward transform 
and thresholding will be zero In the area of the uniform background. 

The convolutional Inverse will therefore receive no contribution from 
outside the retina. As a result one can use the convolutional form of 

the Inverse provided the Image of the ecene le totally contained within 

the ratine. 

2.1.4 Normal Ization. 

Once again one flnde that the reconstructed reflectance le not 
unlqua. That le. any non-singular solution of Llpfa.y)) - 8 can be 
added to tha Input without affecting the result. On the Infinite plane 
auch eolutions have the form p(x t y) » (a*x + b*y + c). If the scene 
only occupies a finite region of space it can be further shown that the 
eolution will be unique up to a constant and that one doee not have to 


worry about possible slopes. To be specifics the background around the 
ecene Ml I I be constant in the reconstruction. So one has here exactly 
the same normalization problem as In the ons-dleenslonal case. 
Assigning uhlts to the region with highest numerical value In the 
reconstructed output appears to be e reasonable method. 

2*1.6 Tuo-Dimeneional Continous (lethodi Details* 

Let r'(x.y) be the reflectance of the object at the point 
corresponding to the iaags point <x t y). Let e* (x,y) be the source 
intensity at that object point. Let pMx.y) be their product, that ie 
the intensity at tht iaage point U.y). Note that rMx.y) and a*(x,y) 
are poei tlve. 

p* (x,y) - a* (x t y) ft p 1 (x,y) 

Let p(x a y) be the logarithm of p* (x.y) and so on. 

p(*.y> - a(x,y) t r(x«y) 

Now assume that s(x,y) and Its first partial derivatives are continoue - 
a reasonable assumption to make for the distribution of illumination on 
the object. Let L be the Laplacian opsrator* 


d(x,y) - L(p(x,y)) - Llsfx.y)) + L(r(x,y)) 

Nom L(i(k»u)I Mill be flnl te everywhere, whl le L(r (x,y) ) mIII bo zero 
except at each edge eeparating regions, where one nil I flr>d a double 
pulse wall ae deecribed. Now discard all finite partet 

t<x,y) - T(L(p<x.y))> -Ltrlx.yM 

Lot G bo the operator corresponding to convolution by (1/2* ) log,(l/r] 
Note that (G)" 1 - L. 

■U,y) - G(T(L(p(* t y))) - r(x t y) +C 

2.2 Tuo-Dlmenslonel flethod - Oiecrete Case. 

Once again we turn fron a contlnous image to one sampled at 
dlecrata points. First Me wit l have to decide on a teaselatlon of the 
Image plane* 

2.2.1 Teaselatlon of the .nage Plans. 

For regular teseelations the choice la between triangular, aquara 
and hexagonal unit eel la* In much past work on image proceeding, aquar 
teaeelations have been used for the obvious reaeone. Thle particular 
taeaalation of the Image has a number of disadvantages. Each cell hae 


two kindt* of neighbors, four adjoining the sides, four on the corners. 
This results In a number of aaymmstr lee. It makes It difficult to find 
convenient difference schemes approximating the Laplacian operator with 
low error term. 

Triangular unit cells are evsn worse In that they have three kinde 
of neighbors* compounded these drawbacks. Note also that near-circular 
objecte pack tightsst in a pattern with hexagonal cella. For these 
reasons hs mIM use a hexagonal unit cell. It should be kept In Bind 
however that It is easy to develop equivalent results using different 
tesselat Ions. 

2.2.2 Discrats Analogue of tha Laplacian. 

Having decidsd on ths teaselation we nssd now to find a discrete 
analogus of the Laplacian operator. Convolution Mlth a central positive 
value and a rotational ly syuetric negative surround of equal weight le 
one possibility. Aside from a negative scale factor, this Mill approach 
application of the Laplacian in the Hull as ths cell size tends to 

If one were to use complicated surrounds, the trade-offs between 
accuracy and resolution would suggest using a negative surround that 
decreasss rapidly outward. For the sake of simplicity we will choose 
convolution with a central cell of weight 1, surrounded by elx eel le of 
weight -1/6. Thle function is convenient, symmstrlc end has a small 
error ten* It Is equsl to - ( h/4 L + h/64 L 1 ) plus sixth and higher 

FIGURE 4a: A discrete analogue of 
the laplacian operator. 

. -*■ 

Nflwvtti Delta function minus 

this discrete analogue. 


order derivative* IRlchtmeyer 19571 . It should again be pointed out 
that similar results can be developed for different functions. 

2-2.3 Inverse of the Oiscrete Opsrator, 

The forward differencing operator has the form 

«»i ■ p., -J}^h * **\ 

Uhere p,« 1* the logarithm of linage intensity, m ( ; fti'e Might*, uhleh In 
our case are 1/G, and the sun is taken over the six immediate neighbors* 

Ue now have to determine the inverse operation that recovers d-- 
from d|i* One approach ie to try and solve the difference equation of 
the form 

v5j Vi ^ p w " -| i 

Or In matrix form U p • d. Note that U la tparse, hevlng l*e on the 
diagonal and -1/G'e scattered around. For a finite retina with n eeneor 
cells one has to introduce boundary conditions to ensure that one has as 
many equations ae there are unknowns. One then eimply inverte the 


matrix U and getst p - U~ d. 

This is entirely analogous to the solution in the continous case 

for a finite retina. U~ corresponds to the Green's function. Much as 
Green's function has a large "support", that is, is non-zero over a 

large area, so 1*1 is not sparse. This implies that a Jot of 
computation is needed to perform the inverse operation. 

2.2.4 Computational Effort and Simplification. 

Solving the difference equations for a given image by simple Gauss- 
Jordan elimination requires of the order of n /2 arithmetic operations. 
Another approach is to invert U once and for all for a given retina. 
For each image then one needs only about n x arithmetic operations. Note 
that the other operations, such as forward differencing, require only 
about Gftn ar i thmet ic operations. 

Uhat in effect is happening is that each point in the output 
depends on each point in the differenced image. Both have n points, so 
n 1 operations are involved. Not only does one have to do a lot of 
computation, but must also store up the matrix U of size n . This le 
quite prohibitive for even a small retina. 

This latter problem can be avoided if one remembers the 
simplification attendant to the use of an infinite retina in the 
continous case. There ne found that the integral ulth Green's function 
simplified into a convolution. Similarly, if one assumes an infinite 
retina here, one finds that W and its inverse become very regular, Tho 


rows in U are then all the sane and the same is true of U . Each value 
in the output then depends in the same way on the neighboring points In 
the differenced image. One need only store up the dependence of one 
point on its neighbors for this simple convolutional operation. 

The only remaining difficulty is that U 19 now Infinite and one can 
no longer invert it numerically - one has to find an analytical 
expression for the inverse. I have not been able to find this inverse 
exactly. A good first approximation is log,lr D /rl - except for r » 0, 
when one uses 1 + log,(r c K Here r is the distance from the origin and 
p^ is arbitrary. The reaainder left over i-ihen one applies the forward 
difference scheme to this approximation lies between log.ll + r ) and 
log 11 - p" ) t This error term is of the order of r • 

In practice one does not have an infinite retina, but as has been 
explained for the continous case one can use the convolutional method 
described above for a finite retina, provided that the image of the 
scene is wholly contained inside the boundaries of the retina* It is 
possible to find an accurate inverse of this kind valid for a limited 
retinal size by numerical means. 

2.2.5 Two-Dimensional Discrete Method: Details. 

Let r*. be the reflectance at the object point corresponding to the 
image point (l.j). Let »y be the intensity of the incident light at 
this object point. Let &l be the intensity in the image at point 


»?j ■•b-4 

Let p;j be the logarithm of pV and so on. Ltt L* be the operator that 
corresponds to convolution with tha analogua of tha Laplaclan* Lat G* 
be i ta Inverse. 

df i * p 'i "2-?^* p w ,d - u( p» J 

Tha weigthe m- are 1/6 in thte casa t and tha eua la taken over the ah 
Immediate nelghbora. 


Here the sum extends over the whole retina and v-» is the convolut ional 
Inveree found numerically as explalnsd above. 


2.2.G Simplicity of the Invert*. 

The forward transform. Involving only a slmpls subtraction of 
Immediate neighbors, Is clearly a rapid, local operation. The inverse 
on the other hand ie global, since each point in the output depends on 
each point in the differenced image. Computationally this makes the 
inverse slow. The inveree is simple in one sense houeveri The 
difference equations bslng solvsd by ths Inverse have the seme form as 
the equations used for the forward transform end are thus local. The 
problem Ie that the output here feeds back into the system and effecte 
can propagate across ths retina. The apparent global nature of the 
inverse is thus of d rather special kind find, as we will see later, 
gives rise to very simple implementations involving only local 
connect tone. 

2.2.7 Iterative Methods of Solution. 

There are of course other methods for solving large sete of 
equations. The fact that U is sparse and has large diagonal elements* 
euggeets trying something I ike Gauss-Ssldsl iteretion. Each Iteration 
takes about S*n arithmetic operations. For effects to propagate acroee 
the retina one requires at least ^(4ftn* - l)/3 iteratlone. Thle le 
becauee a hexagonal retina of width m has (3*m x + l)/4 cells. The above 
euggeets that one might be able to get away with less than n l arithmetic 
operatione. In practice It Is found that effecte propagate very elouly 


and many more Iterations art needed to stabilize the solution, Ona does 
not hava to atora U t alnca it la aaailg generated aa one goea along. 
Iterative echemee correspond to adding a tims-der ivative to tha 
Poiaaon equation and so turning It into a heat-equation, Aa ona 
continues to iterate the steady-state solution la approached. This 
intuitive modal givea some insight into how the process will converge* 

2.2*8 Convergence of Iterative and Feed-back Schemes. 

It is not immediately clear that iterative achamea of solving the 
difference equations hIJI converge. If they do, they will converge to 
the correct solution. Let 5 be the delta function, that is. one at tha 
origin, zero elsewhere. It can be shown that if tha forward 
convolutional operator la u* the convergence of Iterative echemee 
depends on the behaviour of the error term, ( 6 * w) , aa n bacomaa 
large. Raiaing a convolutional operator to an integer power is 
intended to elgnlfy convolution with itself. 

In our case, m is ons at the origin, ulth elx valuea of -1/6 around 
it. So ( 6 - y) ill 1 1 b« zero at the origin with elx valuaa of 1/6 around 
it. New while ( i- w) Ml 1 1 aluays have a total area of one. It doea 
epread out and I te value tends to zero at every point aa n tende to 
Infinity. So thie iterative scheme convsrgeet similar results could bo 
derived for other negative surrounds. 


2.2*9 Setting the Threshold. 

In the dlecrete case a finite threshold must be selected. As 
before, let g' be the seal lest stsp In the logarithm of reflectance In 
the ecene, h the sensor spacing and u the radius of the point-spread 
function of the optical system* Than us define the effective minimum 
step asi 

g • g' * mind, h/2u) 

There are some minor differences In what follows depending on whether 
one considers the sensor outputs to bo intensity aanplee at cell-centore 
or averages over the cell area. The smallest output dus to Mn edge ui II 
be about g/6. This Is producsd uhsn ths edge le oriented to cover Just 
one cell of the neighborhood of six. Let Q be the maximum of the 
intensity gradient - that is the Laplacian of intensity In this case. 
Choose the threshold s as follows: 

e < g/G 

e > * h* 

e > 3^7/G" * 


2*2*10 Sons Notes on This tlethod* 

Notics that an illumination gradient that varies aa soms power of 
distance across ths image becomes a linear slops after taking logarithms 
and thus produces no component after the differencing operations* Such 
simple gradients are eupprseeed even without the thresholding operation. 

In practice ths parameters used in choosing the threshold may not 
be known or may be variable. In this cass ons can look at a histogram 
of the differenced image. It will contain values both positive and 
negative corresponding to edges and also a large number of values 
clustered around zsro due to i I lunination gradients, noise and so on* 
The threshold can bs convsnisntlu choosen to contain this central blob* 
Noise and illumination gradients havs an effect similar to that In 
the one-dimensional cass. Ulth finite cell spacing ons cannot precisely 
eeparate ths two components of ths image intensity and at each edge the 
information Mill be corrupted slightly by nolss and Illumination 
gradient. As the density of sdgas per csll arsa goss up the effect of 
this becomes mors apparent. In highly textured scenes the illumination 
gradient is hard to eliminate* 

Once again one hee to dsclde on a normalization echeme. The best 

method probably Is to 1st the highest numerical value In the 


reconstructed output correspond to unite* 


2.2.11 Dynamic Range Reduction. 

Applying the retinex operation to an image coneiderable reducee the 
range of vatuee. Th'io ie because the output, being related to 
reflectance, will only have a range of one to two ordere of magnitude, 
while the Input will a I eo have i lluminat ion gradients. This will make 
euch processing useful for picture recording and transmission (Horn 

2.2.12 A Frequency Domain Interpretation* 

It may be of Intoreftt to look at thit as t hod fron got another point 
of view. Uhat one does Is to exentuate the high-frequency components, 
threshold and then attenuate the high-frequency coeponente. To eee 
this, coneider first the forward operation. The fourier traneform of 
the convolut tonal operator corresponding to differentiation le 1«. 
Similarly the two-dimensional fourier traneform of the convolutional 
operator corresponding to the LapJeclan le - p*. Here p ie the radiue 
in a polar coordinate system of the two-dimensional frequency space. In 
either case one is multiplying the fourier transform by sone function 
th^t increaeee with frequency. Now consider the revsree operation. The 
fourier traneform of the convolutional operation corresponding to 
Integration ie l/l«. Similarly the fourier transform of (1/2*) log r 
(1/r) ie -l/p x . So In the Inverse step one undoes exactly the emphaele 
given to high frequency components in the forward operation. 


In both thn ons-dlmenslonal and the tuo-dinensional case one loeee 
the zero frequency component. Thle It why the reeult hae to be 
normal ized. 

2.3 Phgelcal llodele. 

There are nueerous continous physical nodele to illustrate the 
inverse traneformatlon. Anything that satisf lee Polsson'a equation ulll 
do. Such physical models help one visualize Mhat the Inverse of a given 
function might be. Examples in two dimensions arei perfect fluld-f Iom» 
eteady diffusion, steady heat-flou, deformation of an elastic membrane, 
electro-statics and current flow in a resistive sheet. In the last 
model for example, the input is the distribution of current flowing Into 
the reeietlve sheet noraal to Its surface, the output ie the 
dletrlbutlon of electrical potential over the surface. 

In addition to helping one visualize solutions* these contlnoue 
eodels aleo suggest discrete models. These can be arrived at elaply by 
cutting up the two-dimensional space in a pattern corresponding to the 
interconnection of neighboring cells. That Is, the remaining parte form 
a pattern dual to that of the sensor cell pettern. Us ulll discuss only 
one such discrete mods It 


2.3.1 A Discrete Physical node I. 

Consider the reeletive eheet described* cut up in the dual pattern 
of the hexagonal unit cell pattern. What will be left ie an 
Interconnection of resistors In a triangular pattern. The inpute to 
thie system will be currents injected at the nodes, the potential at the 
nodes being the output. Thie then provides a very simple analog 
implementation of the tedious inverse computation. 

It ie perhaps at first surprising to sss that each cell ie not 
connected to every other In a direct fashion. One would expect thie 
from the form of the computational Inverse. Each cell In the output 
does Of COUrft* h«V* ft connection Via the ftthtr ceils to oich of thi 
Inputs. Paths are shared however in a way that makes the result both 
simple and planar. 

Conelder for the moment just one node. The potential at the node 
ie the average of the potential of the six nodes connected to it plue 
the current injected times H/B, where R Is the resistance of each 
resistor. The economy of connection is due to the fact that the 
outpute of thie system are fed back into it. It also illuetrates that 
thie model locally solves exactly the same difference equation ae that 
used In the forward transform, only now in reverse. 

This immediately suggests an important property of thie model! By 
eimply changing the Interconnections one can make an inverse for other 
forward transforms. Slmpleet of all are other image plane tesse lat ions, 
both regular and irregular. One simply connecte the resistors in the 

■ '* 

FIGURE 5: Resistive model of the inverse computation. 

The inputs are the currents injected at the nodes. 
The outputs are the potentials at the nodes* 



sane pattern as are tha calls in tha input* 

Mora complicated weighted surrounds can ba handled bg using 
resistors with resistances inversely proportional to tha weights* Tha 
network of resistors will then no longer be planar. 

2*3.2 A Feed-back Scheme for the Inveree. 

Both the comment about outputs feeding back into the resistive 
Model and the earlier notes about iterative schemes suggest get another 
interesting model for the inverss using linear cunning devices. 
Operational amplifiers can serve this purpose. One simply connects the 
summing element so that they solve the difference equation implied by 
the forward transform. Once again It la clear that such a scheme can be 
generalized to arbitrary tesselations and wsighted negative eurrounde 
eimply by changing the Interconnections and attenuations on each input. 
Some questions of stability arise with esotsrlc Interconnections. For 
the simple ones stability Is assured. 

A little thought will show that ths resistive model described 
earlier ie in fact a more economical implementation of just this scheme 
with the difference that there the inputs are currents, while here they 
are potentials. 


























































2.4 Limitations of the Simple Scheme Presented. 

Tha Method presented hara ul II not correctly calculate reflectance 
if uaed unmodified on general scenes. It may however calculate 
lightness fairly well. Ae the msthod stands now for example, a sharp 
ehadow edge ulll not be distinguished from a real edge In the scene and 
the two regions so formed nil 1 produce different outputs, whi le their 
reflectances are the same. It nay be that thie le reasonable 
nevertheleee, elnce ue perceive a difference in apparent lightness. 

Smooth gradations of reflectance on a surface due either to shading 
or variations In surface reflectance ulll be eliminated by the 
thresholding operations except as far as they affect the intensity at 
tha borders of the region* Thie may imply that ue need additional 
channele In our visual eystem to complemsnt the onee cerrylng the 
retlnened Information elnce ue do utilize shading as a depth-cue* 

The simple normalization scheae described Mill also be sensitive to 
specular reflections, flourescent paints and light-sources in the field 
of view* Large depth-discontinul ties present another problem. One 
cannot assume that the illumination is equel on both sides of the 
obscuring edge. In this case tha illuminating component does not vary 
smoothly over the retina, having inetead some sharp edges. 

S -2 

C ** *-> 

3 £ E 



E ^ 

0> r- 

5 2 

'2 i 

E 5 
















2.S Computer Simulation of the Discrete Method* 

A computer program was used to simulate the retinex process 
described on a small retina with both artificial and real images eeen 
through an image dissector camera. The hexagonal unit cell le used in 
this program and the retina itself Is also hexagonal. The retina 
containe 1027 cells in a pattern 36 celle across. Thie ie a compromise 
dictated by the need to limit the number of arithmetic operatlone In tha 
inveree traneform. In this case one needs about a million and thie 
takee about a minute of central processor time on our PDP-18. 

Both the art 1 f Iclal and the real Hondrians consist of regions 
bounded by curved outline* to smphaaiza that this ntthod does not 
require straight-line edges or boundsry extraction and description. 
Various distributions of Incident illumination can be selected for the 
artificial ecenes. In each case the processing satisfactorily removes 
the gradient. 

For the real scenes It is hard to produce really large Illumination 
gradiente by poeltloning the light-sources. The reconstruction does 
eliminate the gradient well, but often minor flaws m1 1 1 appear In the 
output due to noiee in the input and a number of probleme ulth this kind 
of Input device such as a very considerable scatter. It ie not easy to 
predict uhat effects euch imaging device defects will have. 

The output is displayed on a OEC 349 display which has a mere eight 
grey-levele. It uould be interesting to experiment ulth larger retinae 
and better Image Input- and output-devices. 


2.5.1 For« of Inverse used In the Computer Simulation* 

The convo I u t i one I form of the Inverse was used for speed and lou 
storage requirement. This neccesitated solving the difference equation! 
once, given a pulse as Input. The symmetry of the hexagonal pattern 
allous one to identify symmetrically placed cells and only 324 unknowns 
needed to be found for a convolutional inverse sufficient for the eize 
of rotina described. As Mentioned before, this function is closely 
approximated by log, (r # /rl for large r. This can be used to establish 
boundary condi t ions. 

FIGURE 8: The method applied to an 
artificial image. 

FIGURE 9: The method applied to a 
real image 

FIGURE 10: The method applied to 
Craik's figure. 

FIGURE 11: Apparent lightness predicted 
for incomplete figure- 

The subfigures in the above have the following interpretation: 

A Input - logarithm of image intensity 

B Differenced image 

C Thresholded difference 

D Output - computed lightness 

E Illumination distribution (P^*!**) 





3* Implications and Conclusions. 
3*1 Parallel Imaga Procsaaing Hardware. 

Tha methods described hare for forward transforming, thresholding 
and Inverse transforming inaediately tempt one to think In terms of 
electronic componenta arranged in parallel layere. Enough has been aald 
about different models to make it clear how one Bight connect euch 
components. Large scale integrated circuit technology may be useful* 
provided the signals are either convsrtsd from analog to digital for* or 
batter still, good linear circuits are available In thia form. 

Construction of such devlcea would be premature until further 
experimentation le performed to decide on optimal tesselatlone, optimal 
negative surrounds, thresholding operations and normalization schemes. 
Theee decleions are best guided bg computer simulation. 

3.2 Cognitive Peychology. 

One of the artificial ecanes was created to Illustrate Cralk's 
lllueion (Brindley I960, Cornsweet 1970). Here a sharp edge le bordered 
by eecond-order gradients. As one might expect the smooth gradlente are 
loet in the thresholding and reconstruction produces two regions each of 
uniform brightness* Ths diffsrencs in brightnsss between the reglone le 
equal to the original Intensity stsp at the edge. 

The fact that the process presented hsre felle prey to this 


i I lueion is of course no proof that humans use the «a»* mechanism. It 
le interesting that this technique allows one to predict for example the 
appearance of pictures containing incompletely closed curves with 
second-order gradients on either side. 

3*3 Neuro-physr ology. 

The method described here for obtaining lightness from Image 
Intmnelty euggeete functions for a number of etructuree in the primate 

retina. The horizontal ceils appear to be involved In the forward 
transformation, while some of the amacrlnes may be Involved in the 
inverse transformation* For details see tha paper by David Narr IDarr 
19741 , In which he uses this hypothssis to explain an astonishing number 

of facte about the retina. 

3.4 Concluelon. 

A elmple layered, parallel technique for computing lightness from 
image Inteneity has been preesntsd. Ths method does not involve an 
ability to describe or understand the ecene, relying instead on the 
spatial differences In the distribution of reflectance and illumination. 
The forward step involves accentuating ths edges between regions. The 
output of thle step is then thresholded to remove Illumination gradiente 
and noise. The Inverse step merely undoes the accentuation of the 


Physical models nave been glvan uhich can perform this commutation 
efficiently in parallel layers of simple networks. The method haa been 
simulated and applied to a number of images* The method grew out of an 
attempt to extend Land'e method to two dimenslone and fills the need for 
a llghtneee Judging process in his retinex theory of color perception. 

The possibility of processing an imags In such a parallel* elmpla 
fashion uithout higher-level understanding of the scene reinforces my 
belief that such lou-level processing is of Importance In dealing ulth a 
number of features of Images. Amongst thsee are shading, stereo 
disparity, focue, edge detection, scene segmentation and motion 
parallax. Some of thia kind of processing mag actually happen In the 
primate retina and viaual cortex. The Implications for Image analysis 
are that it may well be that a number of such pre-proceeeing operatlona 
should be performed automatically for the whols Image to accentuate or 
extract certain attributes before one brings to bear the more powerful* 
but tedious and slou sequential goal-directed methods. 


B lb I iographg. 

Brindleu. G.S, (1368) "Physiology of the Retina and Visual Pathway." 
(Honograph No. 6 of the Physiological Society). Londoni Edward 
Arnold Ltd. 

Cornsueet, T. (1978) "Visual Perception." Nau Vorki Academic Preae. 

Garabedian, O.R. (19G4) 'Partial Olffarantlal Equations." Nau Yorki John 
Ui lay. 

Goetha, J. 14. von. (1818) "Zur Farbenlehre. " Tuebingen. 

Hardy, A.C. (Ed,) (193S) "Tha Handbook of Coioriaetry. " Cambridge, Haest 
n. I ,T» Praaa* 

Hefmhoitz, H.L.F. (1867) "Handbuch dar Physlologlachan Optlk." Lalpzlgi 
Voaa. Alao tranalatedi Southall, J.P.C.S. "Handbook of 
Physiological Optics." Nau Vorki Oovar Publications. 

He I son, H. (1938) "Fundamental Problems In Color Vision I." Journal of 
Experimental Psychology, 23. 

Helaon, H. (1948) "Fundamental Problems in Color Vision II." Journal of 
Experimental Psychology, 26. 

Herlng, E. (1875) "Zur Lahrs von Lichtsinne - Grunzuege elner Theorie 
des Farbsinnes." SBK. Akad. Uiss. Uien ftath. Naturuies. K.7B. 

Horn, B.K.P, (1368) "Tha Application of Fourier Transform Methods to 

Image Procasalng." S.M. Thesis. E.E. Department, fl.I.T. pp. 81-87 A 
pp. 93-94. 

Judd, D.B* (1948) "Hue, Saturation and Lightness of Surface Colore with 
Chromatic Illumination ." Journal of the Optical Society of 
America, 38. 

Judd, D.B. (1962) "Color In Business, Science and Industry." New Vorki 
John Ulley. 

Land, E.H. (1959) "Experlmsnte in Color Vieion," Scientific American. 

Land, E.H, (1964) "The Retinex," American Scientist. S2. 

Land, E.H. 4 flcCann, J.J. (1971) "Lightness Theory." Journal of the 
Optical Society, 61. 

Lettvin, J.Y. (1967) 'The Color of Colored Things." Quarterly Progress 
Report, 87, Rssearch Laboratory for Electronics, M.I.T. 


Marr f D- (1974) "An Analysis of the Primate Retina." A.I. Clemo, 296, 
Artificial Intel I Igence Laboratory, M.I.T. 

tlaxueH, J.C. (185G) "On the Unequal Sensitivity of the Foramen Centrale 
to Light of Different Colours. " Rep. Brit. Assoc. 

Newton, Sir Isaac. (1784) "OpticKs." Londoni Sanuel Smith 6 Benjamin 
Uatford. Also Neu Yorkt Dover Publ icatione. 

Richards, U. (1971) "One-stage flodsl for Color Convereion. ■ Journal of 
the Optical Society, 62. 

Rlchtmeyer, R.O. & dor ton, K.U. (1957) "Difference Hethode for Initial 
Value Problems.- Neu YorKi John U isy. 

Young, T. (1828) "On the Theory of Light and Colour." Philosophical 

Transact lone* Aleo im Teevan, R.C. & Birney, R.C. (Ed.) (1361) 
"Color Vision," Van Noetrand.